TempRanges

Code

spec_df<-data.frame(sci=c("Epinephelus morio","Lutjanus campechanus","Mycteroperca microlepis","Mycteroperca phenax","Rachycentron canadum","Rhomboplites aurorubens","Seriola dumerili"),
                    com=c("Red Grouper", "Red Snapper", "Gag Grouper","Scamp Grouper","Cobia","Vermilion Snapper","Greater Amberjack"))

scinm<-spec_df$sci
comnm<-spec_df$com

DATA

GOAL: We wanted a better way to examine species temperature ranges based off of real data, take ranges from a distribution to combat sampling bias, a way to semi “ground truth” general distributions, and a way to inform the depth of the area of analysis. The data we have available to us are the video surveys with related temps (bottom temp of survey?).

Tip

You can click on figures in this document to see them enlarged and zoom in!

FYI: For those interested this is possible by using the lightbox Quarto extension.

Code

#2011-2017
mosttemps<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/AllSpeciesSumcounts2011-2017.csv")
mosttemps<-mosttemps %>% dplyr::select(Year, Collection, LastOfTemp)

#2018
temps18<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/videopresenceabsence2015-2018.csv")
temps18<-temps18 %>% dplyr::select(Year, Collection, LastOfTemp)
temps18<-temps18[temps18$Year>2017,]

#All
temps<-rbind(mosttemps, temps18)
temps<-na.omit(temps)

Code

vid_ab<-read_xlsx("C:/Users/brittanytroast/Dropbox/Work/Random/Before Fixing Comp/Climate Fisheries Initiative/Data/Video_Counts_2011-2018.xlsx")

vid_ab<-vid_ab[,c(1,2,4,47:ncol(vid_ab))]
vidab_long<-pivot_longer(vid_ab, c(4:ncol(vid_ab)), names_to = "Species", values_to = "AB")

vid_spec_ab<-vidab_long[vidab_long$Species %in% scinm,]
# unique(vid_spec_abna$Species)
vid_spec_abna<-vid_spec_ab[vid_spec_ab$AB>0,]

Code

colnames(temps)[3]<-"Temp"
temps$Year<-NULL
vid_spec_abna$Collection<-as.integer(vid_spec_abna$Collection)
vid_abtemp<-left_join(vid_spec_abna, temps, by=c("Collection"))

vid_abtemp<-na.omit(vid_abtemp)
vid_abtemp<-merge(vid_abtemp, spec_df, by.x="Species", by.y="sci")

TEMP

I first “squinted” at the temperature data per species to determine general curves. Here are the temperature distributions for the sum of catch at half degree temperature steps per species.

Code

#manual binning SUM
vid_bin<-vid_abtemp %>% mutate(tempbin=cut(Temp, breaks = seq(16, 31,0.5)), tempmid=as.numeric(str_sub(gsub("^(.*?),.*", "\\1", tempbin), 2))+.25)
vid_sum<-vid_bin %>% group_by(Species, tempmid) %>% summarise(sum=sum(AB))

vid_sum<-merge(vid_sum, spec_df, by.x="Species", by.y="sci")
ggplot(vid_sum, aes(x=tempmid, y=sum))+
  facet_wrap(~com, scales = "free_y", ncol=2)+
  # geom_bar(stat="identity", aes(fill=Species), color="black")+
  stat_summary(geom="bar", fun="sum", aes(fill=Species), color="black")+
  geom_smooth( method="loess", color="black")+
  scale_y_continuous(limits = c(0,NA))+
  labs(title = "Sum Abundance 0.5°C 2015-2017 SA",x="Temperature °C", y="Sum CPUE")+
  theme_bw()+theme(legend.position = "none",
                   strip.background = element_blank(),
                   strip.text = element_text(face="bold"))

Code

ggsave("AB_HalfDeg_Sum.png", width = 6, height = 6, units="in")

To address the data richness bias here are the 10-90% percentile temp ranges. This table includes these values as well as min and max temps the species was sampled at. This data is from the video trap survey data from 2011-2018 in the U.S. South Atlantic, South of Cape Hatteras, North Carolina. While it is not the whole area we are interested in, it does at least give us an idea of temperature ranges. If we had the northern data we might be able to get better cold temperature data, but luckily we care more about the upper temperature limit.

Code

sa_temp<-vid_abtemp %>% group_by(Species, com) %>% summarise(mintemp=round(min(Temp),2),ten_quant=round(quantile(Temp,prob = c(.10)),2),median=median(Temp), ninquant=round(quantile(Temp,prob = c(.90)),2),maxtemp=round(max(Temp),2))

# sa_temp %>% gt(rowname_col = "com", groupname_col = NULL) %>%
#   tab_header(title="Temperature Range °C") %>%
#   cols_label(Species="Scientific Name", mintemp="Lowest Temp",ten_quant="Lower Quantile", ninquant="Upper Quantile", maxtemp="Highest Temp") %>%
#   tab_stubhead(label="Common Name") %>%
#   tab_style(style=cell_text(style = "italic"), locations = cells_body(columns = Species)) %>%
#   tab_style(style = cell_text(color = "darkgrey", 
#         font = google_font("Source Sans Pro"), transform = "uppercase"), 
#         locations = cells_stubhead()) %>%
#   gt_theme_nytimes()


sa_temp<-sa_temp %>% mutate(comb=paste0("<b>",com, "</b><br><i>", Species,"</i>"), ten_quant=round(ten_quant, 2), ninquant=round(ninquant, 2))


sa_temp %>% gt(rowname_col = "comb", groupname_col = NULL) %>%
  cols_hide(columns=c(Species, com, median)) %>%
  tab_header(title="Temperature Range °C") %>%
  cols_label(mintemp="Lowest Temp",ten_quant="Lower Quantile", ninquant="Upper Quantile", maxtemp="Highest Temp") %>%
  tab_stubhead(label="Common Name") %>%
  tab_style(style=cell_text(style = "italic"), locations = cells_body(columns = comb)) %>%
  tab_style(style = cell_text(color = "darkgrey", 
        font = google_font("Source Sans Pro"), transform = "uppercase"), 
        locations = cells_stubhead()) %>%
  gt_theme_nytimes() %>%
  fmt_markdown()

Temperature Range °C
Common Name	Lowest Temp	Lower Quantile	Upper Quantile	Highest Temp
Red Grouper Epinephelus morio	15.16	18.51	25.71	28.06
Red Snapper Lutjanus campechanus	12.92	18.76	26.08	29.31
Gag Grouper Mycteroperca microlepis	13.62	19.11	26	28.25
Scamp Grouper Mycteroperca phenax	13.05	18.55	26.22	28.4
Cobia Rachycentron canadum	13.56	20.02	25.63	28.66
Vermilion Snapper Rhomboplites aurorubens	13.05	18.95	26.7	29.31
Greater Amberjack Seriola dumerili	13.71	19.63	26.64	29.26

I further explored the data to see how many data points (each point represents one individual caught at a give temperature) we were excluding using the 10-90% range. Unlabeled text are max temperature, 90% percentile temperature, 10% percentile temperature and minimum temperature respectively from top to bottom. Labels within the boxplot show the total individuals caught (Total n), samples where individuals occurred (Sample n), and the number of individuals excluded by the 10-90% range (Out n). I figured we could use this to help us refine which range we would like to use (e.g. if we were excluding large numbers of individuals maybe we should expand our range or find natural breaks).

Code

#Single instances of abundance (for jitter)
vid_exp<-vid_abtemp %>% uncount(AB) %>% mutate(count=1)

#Add points outside of 10-90% to show how many are outside that range
temp_li<-list()
specs<-spec_df$com
for (i in 1:length(spec_df$com)) {
  dat<-vid_exp[vid_exp$com==specs[i],]
  lowup<-sa_temp[sa_temp$com==specs[i],]
  dat$inout<-ifelse(dat$Temp < lowup$ten_quant | dat$Temp > lowup$ninquant, "out", "in")
  temp_li[[i]]<-dat
 }
inout_df<-do.call("rbind", temp_li)

#Add Total n & Sampling n
ndf_1<-vid_abtemp %>% group_by(Species, com) %>% summarise(tot_n=sum(AB), samp_n=length(AB))
ndf_2<-inout_df[inout_df$inout=="out",] %>% group_by(Species, com) %>% summarise(out_n=sum(count))
n_df<-left_join(ndf_1, ndf_2)
n_df$lab<-paste0("Total n=", n_df$tot_n, "\nSample n=", n_df$samp_n, "\nOut n=", n_df$out_n)


#Boxplots for data
ggplot(sa_temp, aes(x=com))+
  geom_jitter(data=inout_df[inout_df$inout=="out",], aes(y=Temp),color="gray60", alpha=0.3, width = 0.4,height = 0.2)+
  geom_boxplot(stat="identity",aes(ymin=mintemp, lower=ten_quant, middle=median,upper=ninquant, ymax=maxtemp, fill=com))+
  geom_text(aes(label=mintemp, y=mintemp-0.5), size=3.5, fontface="bold")+
  geom_text(aes(label=maxtemp, y=maxtemp+0.5), size=3.5, fontface="bold")+
  geom_text(aes(label=ten_quant, y=ten_quant+0.5),color="white", fontface="bold",size=3.5)+
  geom_text(aes(label=ninquant, y=ninquant-0.5),color="white", fontface="bold", size=3.5)+
  geom_label(data=n_df, aes(label=lab, y=sa_temp$median), size=2.5, fontface="bold")+
  labs(title = "Video Survey Temp 10-90% Quantile" , x="Species", y="Temperature")+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+theme(legend.position = "none",
                   title = element_text(face="bold"))

Unlabeled text are max temperature, 90% percentile temperature, 10% percentile temperature and minimum temperature respectively from top to bottom. Labels within the boxplot show the total individuals caught (Total n), samples where individuals occurred (Sample n), and the number of individuals excluded by the 10-90% range (Out n).

Based on the above figure I explored a broader range of 5-95% and see if it was appropriate.

Code

#| fig-cap: Unlabeled text are max temperature, 90% percentile temperature, 10% percentile temperature and minimum temperature respectively from top to bottom. Labels within the boxplot show the total individuals caught (Total n), samples where individuals occurred (Sample n), and the number of individuals excluded by the 10-90% range (Out n). 

low_quant=0.05
high_quant=0.95

sa_temp_alt<-vid_abtemp %>% group_by(Species, com) %>% summarise(mintemp=round(min(Temp),2),ten_quant=round(quantile(Temp,prob = c(low_quant)),2),median=median(Temp), ninquant=round(quantile(Temp,prob = c(high_quant)),2),maxtemp=round(max(Temp),2))


#Single instances of abundance (for jitter)
vid_exp<-vid_abtemp %>% uncount(AB) %>% mutate(count=1)

#Add points outside of 10-90% to show how many are outside that range
temp_li<-list()
specs<-spec_df$com
for (i in 1:length(spec_df$com)) {
  dat<-vid_exp[vid_exp$com==specs[i],]
  lowup<-sa_temp_alt[sa_temp_alt$com==specs[i],]
  dat$inout<-ifelse(dat$Temp < lowup$ten_quant | dat$Temp > lowup$ninquant, "out", "in")
  temp_li[[i]]<-dat
 }
inout_df<-do.call("rbind", temp_li)

#Add Total n & Sampling n
ndf_1<-vid_abtemp %>% group_by(Species, com) %>% summarise(tot_n=sum(AB), samp_n=length(AB))
ndf_2<-inout_df[inout_df$inout=="out",] %>% group_by(Species, com) %>% summarise(out_n=sum(count))
n_df<-left_join(ndf_1, ndf_2)
n_df$lab<-paste0("Total n=", n_df$tot_n, "\nSample n=", n_df$samp_n, "\nOut n=", n_df$out_n)


#Boxplots for data
ggplot(sa_temp_alt, aes(x=com))+
  geom_jitter(data=inout_df[inout_df$inout=="out",], aes(y=Temp),color="gray60", alpha=0.3, width = 0.4,height = 0.2)+
  geom_boxplot(stat="identity",aes(ymin=mintemp, lower=ten_quant, middle=median,upper=ninquant, ymax=maxtemp, fill=com))+
  geom_text(aes(label=mintemp, y=mintemp-0.5), size=3.5, fontface="bold")+
  geom_text(aes(label=maxtemp, y=maxtemp+0.5), size=3.5, fontface="bold")+
  geom_text(aes(label=ten_quant, y=ten_quant+0.5),color="white", fontface="bold",size=3.5)+
  geom_text(aes(label=ninquant, y=ninquant-0.5),color="white", fontface="bold", size=3.5)+
  geom_label(data=n_df, aes(label=lab, y=sa_temp_alt$median), size=2.5, fontface="bold")+
  labs(title = paste0("Video Survey Temp ",low_quant*100,"-",high_quant*100,"% Quantile") , x="Species", y="Temperature")+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+theme(legend.position = "none",
                   title = element_text(face="bold"))

Lastly, I compared potential new ranges with what we had originally from a lit search.

Code

#Add old upper/lower temps to spec_df
spec_df$lowtemp=c(16,19,18.3,20.7,16.8,18.3,16.9)
spec_df$hightemp=c(25.7,29,27.1,27.5,32,27.2,29)

#Create df with old, 10-90%, and 5-95% quantile temps
all_temps<-data.frame(spec=spec_df$sci, com=spec_df$com,
                      low_old=spec_df$lowtemp, upp_old=spec_df$hightemp,
                      low_1090=sa_temp$ten_quant, upp_1090=sa_temp$ninquant,
                      low_595=sa_temp_alt$ten_quant, upp_595=sa_temp_alt$ninquant)

all_temps<-all_temps %>% mutate(comb=paste0("<b>",com, "</b><br><i>", spec,"</i>"))


all_temps %>% gt(rowname_col = "comb", groupname_col = NA) %>%
  cols_hide(columns=c(spec, com)) %>%
  tab_header(title="Temperature Range Comparisons") %>%
  tab_spanner(.,label = md("**Old**"), columns=ends_with("_old")) %>%
  tab_spanner(.,label = md("**10-90%**"), columns=ends_with("_1090")) %>%
  tab_spanner(.,label = md("**5-95%**"), columns=ends_with("_595")) %>%
  cols_label(starts_with("low")~"Lower", starts_with("upp")~"Upper") %>%
  tab_stubhead(label="Species") %>%
  tab_style(style=cell_text(style = "italic"), locations = cells_body(columns = comb)) %>%
  tab_style(style = cell_text(color = "darkgrey", 
        font = google_font("Source Sans Pro"), transform = "uppercase"), 
        locations = list(cells_stubhead(), cells_column_labels())) %>%
  tab_style(style = cell_text(align = "center", weight = "bold"),locations = cells_title()) %>%
  tab_style(
    style = cell_borders(
      sides = c("right"),
      weight = px(2),
      color = "#d3d3d3"),
    locations = cells_body(
      columns = c(upp_old, upp_1090))) %>%
  gt_theme_nytimes() %>%
  fmt_markdown()

Temperature Range Comparisons
Species	Old		10-90%		5-95%
	Lower	Upper	Lower	Upper	Lower	Upper
Red Grouper Epinephelus morio	16	25.7	18.51	25.71	16.82	26.95
Red Snapper Lutjanus campechanus	19	29	18.76	26.08	17.58	27.18
Gag Grouper Mycteroperca microlepis	18.3	27.1	19.11	26	17.24	26.71
Scamp Grouper Mycteroperca phenax	20.7	27.5	18.55	26.22	17.37	27
Cobia Rachycentron canadum	16.8	32	20.02	25.63	19.35	26.63
Vermilion Snapper Rhomboplites aurorubens	18.3	27.2	18.95	26.7	17.57	27.41
Greater Amberjack Seriola dumerili	16.9	29	19.63	26.64	18.59	27.2

Code

#Plot to visualize temp ranges easily
all_temps<-all_temps %>% pivot_longer(cols = -c(spec, com, comb)) %>% separate(name, c("temp_cat","sceno"))
all_temps$sceno<-factor(all_temps$sceno, levels = c("old","1090","595"))
temp_range<-all_temps %>% pivot_wider(names_from = temp_cat, values_from = value)
temp_range$com_nm<-gsub( " ", "\\\n",temp_range$com)
temp_range$com_nm<-factor(temp_range$com_nm,levels=c("Red\nGrouper","Red\nSnapper","Gag\nGrouper",'Scamp\nGrouper',"Cobia","Vermilion\nSnapper","Greater\nAmberjack"))


ggplot(temp_range, aes(x=com_nm, group=sceno))+
  geom_linerange(aes(ymin=low, ymax=upp, color=sceno, lty=sceno), position = position_dodge(width = -0.6), size=1.25, alpha=0.5)+
  geom_point(aes(y=low, color=sceno), position = position_dodge(width = -0.6), size=3, shape=21, fill="white", stroke=0.9)+
  geom_point(aes(y=upp, color=sceno), position = position_dodge(width = -0.6), size=3, shape=21, fill="white", stroke=0.9)+
  scale_color_manual(values = c("black", "red", "blue"), labels=c("Old", "10-90%", "5-95%"))+
  scale_linetype_manual(values = c("solid", "longdash","dashed"),labels=c("Old", "10-90%", "5-95%"))+
  labs(y="Temperature", x="Species", color="Distribution", lty="Distribution")+
  scale_y_continuous(breaks = seq(16,32, 2))+
  coord_flip()+
  scale_x_discrete(limits=rev)+
  theme_bw() + theme(axis.text.y = element_text(face="bold"),
                     legend.position = c(0.933,0.1),
                     legend.background = element_rect(color="black"))

Code

ggsave("CompareTemps.png", width=8, height=6, units="in")

DISTRIBUTION

Here I explore the general distribution based on the video trap survey in the southeast. I’m not sure this is the best data we could use but again it gives a general idea.

Code

#2011-2017
mostcoords<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/AllSpeciesSumcounts2011-2017.csv")
mostcoords<-mostcoords %>% dplyr::select(Year, Collection, Start_Latitude, Start_Longitude)

#2018
coords18<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/videopresenceabsence2015-2018.csv")
coords18<-coords18 %>% dplyr::select(Year, Collection, Start_Latitude, Start_Longitude)
coords18<-coords18[coords18$Year>2017,]

#All
coords<-rbind(mostcoords, coords18)
coords<-na.omit(coords)

This figure gives a smoothed heat map of density of catches. The data here are aggregated from 2011-2018. This figure is easier to use than the next one but can smooth over and hide some information.

Code

#Raw Data
vidab<-left_join(vid_spec_abna, coords)
vid_ab<-dplyr::select(vidab, year=Year, spec=Species, count=AB, lat=Start_Latitude, long=Start_Longitude)
vid_ab<-left_join(vid_ab, spec_df, by=c("spec"="sci"))

#Coastline
atl<-readOGR("C:/Users/brittanytroast/Dropbox/Work/Random/Before Fixing Comp/Climate Fisheries Initiative/Stats/CFI Hab Envelopes/Shapefiles/AtlCoastBuff.shx", verbose = F)
atl_sf<-st_as_sf(atl)
ext_atl<-extent(atl)
land<-maps::map("usa", fill=T, plot=F)
land_base<-st_as_sf(land)
land_base<-st_make_valid(land_base)
sf_use_s2(FALSE)
land_base_val<-st_make_valid(land_base)
land_atl<-st_crop(land_base_val, xmin =ext_atl@xmin, xmax=ext_atl@xmax, ymin=ext_atl@ymin-2, ymax=ext_atl@ymax)

This figure is similar to the previous however has hex bins with the number of individuals caught from 2011-2018. This smooths over less of the data.

Code

vid_ab$com<-gsub(" ", "\\\n", vid_ab$com)
ggplot(vid_ab, aes(x=long, y=lat))+
  facet_wrap(~com)+
  geom_sf(data=land_atl, inherit.aes = F)+
  stat_density_2d(aes(fill=..level..), geom="polygon")+
  labs(title = "Aggregrate Catch 2011-2018", fill="Level")+
  scale_fill_continuous(type = "viridis") +
  xlim(-82,-75) + ylim(25,36)+
  theme_bw()+ theme(strip.background = element_blank(),
                    strip.text = element_text(face="bold"),
                    title = element_text(face="bold"),
                    axis.title = element_blank(),
                    legend.position = c(0.85,0.15),
                    legend.title.align=0.5)

Code

ggplot(vid_ab, aes(x=long, y=lat))+
  facet_wrap(~com)+
  geom_sf(data=land_atl, inherit.aes = F)+
  geom_hex(bins = 20) +
  labs(title = "Aggregrate Catch 2011-2018", fill="Abundance")+
  scale_fill_continuous(type = "viridis") +
  xlim(-82,-75) + ylim(25,36)+
  theme_bw()+ theme(strip.background = element_blank(),
                    strip.text = element_text(face="bold"),
                    title = element_text(face="bold"),
                    axis.title = element_blank(),
                    legend.position = c(0.85,0.15),
                    legend.title.align=0.5)  

Here is the same heat map as above except with yearly information per species.

Code

dir.create("Spec Heat Maps")
vid_ab<-vid_ab %>% mutate(lat=round(lat, 2), long=round(long,2))

spec<-unique(vid_ab$spec)

for (i in 1:length(spec)) {
  sp_df<-vid_ab[vid_ab$spec==spec[i],]
  sp_df<-na.omit(sp_df)
  p<-ggplot(sp_df, aes(x=long, y=lat))+
  facet_wrap(~year)+
  geom_sf(data=land_atl, inherit.aes = F)+
  stat_density_2d(aes(fill=..level..), geom="polygon")+
  labs(title = unique(sp_df$com), fill="Level")+
  scale_fill_continuous(type = "viridis") +
  xlim(-82,-75) + ylim(25,36)+
  theme_bw()+ theme(strip.background = element_blank(),
                    strip.text = element_text(face="bold"),
                    title = element_text(face="bold"),
                    axis.title = element_blank(),
                    legend.position = c(0.85,0.15),
                    legend.title.align=0.5)
  filename<-paste0(getwd(),"/Spec Heat Maps/HeatMap_",substr(spec[i],1,3),"_",substr(sub(".*? ", "", spec[i]),1,3),".png")
  print(p)
# ggsave(plot = p, filename = filename, width = 5, height = 7, units = "in")

}

Here is the same hex bin map as above except with yearly information per species.

Code

dir.create("Spec Hex Maps")
vid_ab<-vid_ab %>% mutate(lat=round(lat, 2), long=round(long,2))

spec<-unique(vid_ab$spec)

for (i in 1:length(spec)) {
  sp_df<-vid_ab[vid_ab$spec==spec[i],]
  sp_df<-na.omit(sp_df)
  
  p<-ggplot(sp_df, aes(x=long, y=lat))+
  facet_wrap(~year)+
  geom_sf(data=land_atl, inherit.aes = F)+
  geom_hex(bins = 20) +
  labs(title = unique(sp_df$com), fill="Abundance")+
  scale_fill_continuous(type = "viridis") +
  xlim(-82,-75) + ylim(25,36)+
  theme_bw()+ theme(strip.background = element_blank(),
                    strip.text = element_text(face="bold"),
                    title = element_text(face="bold"),
                    axis.title = element_blank(),
                    legend.position = c(0.85,0.15),
                    legend.title.align=0.5)  
  filename<-paste0(getwd(),"/Spec Hex Maps/HexMap_",substr(spec[i],1,3),"_",substr(sub(".*? ", "", spec[i]),1,3),".png")
  print(p)
# ggsave(plot = p, filename = filename, width = 5, height = 7, units = "in")

}

DEPTH

We additionally discussed looking into the depth range of species so that we can customize our area of analysis to better represent an area that species may be able to inhabit if the bottom temperature is suitable.

I did a quick google/lit search of depth limits per species. Some were reputable sources like the fishing council and some were wiki, but I can dive deeper when we determine a game plan.

Code

depth_df<-data.frame(sci=c("Epinephelus morio","Lutjanus campechanus","Mycteroperca microlepis","Mycteroperca phenax","Rachycentron canadum","Rhomboplites aurorubens","Seriola dumerili"),
                    com=c("Red Grouper", "Red Snapper", "Gag Grouper","Scamp Grouper","Cobia","Vermilion Snapper","Greater Amberjack"),
                    depth_upp=c(5,9,18,0,0,18,18),
                    depth_low=c(330,189,76,100,76,122,73))

depth_df<-depth_df %>% mutate(comb=paste0("<b>",com, "</b><br><i>", sci,"</i>"))

depth_df %>% gt(rowname_col = "comb", groupname_col = NULL) %>%
  cols_hide(columns=c(sci, com)) %>%
  tab_header(title="Depth Range(m)") %>%
  cols_label(depth_upp="Upper Depth",depth_low="Lower Depth") %>%
  tab_style(style=cell_text(align="center"), locations = cells_body()) %>%
  gt_theme_nytimes() %>%
  fmt_markdown()

Depth Range(m)
	Upper Depth	Lower Depth
Red Grouper Epinephelus morio	5	330
Red Snapper Lutjanus campechanus	9	189
Gag Grouper Mycteroperca microlepis	18	76
Scamp Grouper Mycteroperca phenax	0	100
Cobia Rachycentron canadum	0	76
Vermilion Snapper Rhomboplites aurorubens	18	122
Greater Amberjack Seriola dumerili	18	73

Here is the above data represented in the same way as the temperature ranges previously so we can compare among species.

Code

depth_df$com<-factor(depth_df$com, levels = unique(depth_df$com))
ggplot(depth_df, aes(x=com))+
  geom_linerange(aes(ymin=depth_upp, ymax=depth_low), size=1.25, alpha=0.5)+
  geom_point(aes(y=depth_upp, fill=com), size=3, shape=21, stroke=0.9)+
  geom_point(aes(y=depth_low, fill=com), size=3, shape=21, stroke=0.9)+
  labs(y="Depth (m)", x="Species", title="Depth range of Selected Species (lit search)")+
  coord_flip()+
  scale_x_discrete(limits=rev)+
  theme_bw() + theme(axis.text.y = element_text(face="bold"),
                     legend.position = "none",
                     legend.background = element_rect(color="black"))

I looked at the depth provided from the video trap survey but examining the only depth metric they have while sampling, it appears they do not sample even close to the full range of some of these species so I don’t think we can use that data.

Code

#2011-2017
mostdepth<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/AllSpeciesSumcounts2011-2017.csv")
mostdepth<-mostdepth %>% dplyr::select(Year, Collection, Start_Depth)

#2018
depth18<-read.csv("C:/Users/brittanytroast/Dropbox/Work/Climate Fisheries Initiative/Data/Full Temps/videopresenceabsence2015-2018.csv")
depth18<-depth18 %>% dplyr::select(Year, Collection, Start_Depth)
depth18<-depth18[depth18$Year>2017,]

#All
depth<-rbind(mostdepth, depth18)
depth<-na.omit(depth)


ggplot(depth, aes(x=Start_Depth))+
  geom_histogram(fill="#3788A6", color="black")+
  labs(title="Depth of Video Trap Sampling")+
  theme_bw()+ theme()

In order to refine the shape of our area of analysis, I tried to find shapefiles of the U.S. continental shelf. For some reason I could not find this easily. I did however find that around 140m is often where the cutoff occurs (shallow continental shelf). So instead I learned how to make my own shapefiles based on whatever depth I would like and found a really nice bathymetry tool/package that uses the ETOPO 2022 database hosted on the NOAA website.

Here is the bathymetric map of the U.S. Atlantic coast. Here I have added contour lines for 140m (general shallow continental shelf), 200m (second deepest range of our species), and 350m (deepest depth limit). I was surprised to find the area difference between these three depths is minimal.

Code

#| fig-cap: Depth is presented in a blue gradient from 0 to 6000 meters. Depth contour lines at 140m (light gray solid line), 200m (medium gray longdashed line), and 350m (dark grey dashed line).

#Get Bathy
bathy<-getNOAA.bathy(lon1 = -82,lon2 = -62, lat1 = 24.9, lat2 = 45)
bath_fort<-fortify.bathy(bathy)
bath_fort<-bath_fort[bath_fort$z<0,]
# bath_fort<-bath_fort[bath_fort$z> -2000,]
#Cut to 140m
# bath_cut<-bath_fort[bath_fort$z<0 & bath_fort$z> -140,]

# blue.col <- colorRampPalette(c("darkblue", "lightblue"))
# library(metR)

ggplot(bath_fort, aes(x=x, y=y, z=z))+
  geom_raster(aes(fill=z))+
  geom_contour(aes(color=c(as.factor(..level..)), lty=c(as.factor(..level..))), breaks = c(-140,-200,-350))+
  geom_sf(data=land_atl, inherit.aes = F, fill="gray70")+
  scale_fill_viridis_c(option = "G")+
  scale_color_manual(values = c("gray25", "gray45","gray65"), labels=c("350m", "200m","140m"))+
  scale_linetype_manual(values= c("dashed", "longdash", "solid"), labels=c("350m", "200m","140m"))+
  guides(fill=guide_colorbar(order=2), color=guide_legend(order=1,keyheight = 2, reverse=T), lty=guide_legend(order=1,keyheight = 2, reverse=T))+
  labs(title="Bathymetry", x=NULL, y=NULL, fill="Depth(m)", color="Contour", lty="Contour")+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  coord_sf(xlim = c(-82,-62), ylim = c(24.9,45))+
  theme_bw()+theme(legend.position = c(0.833,0.134),
                   title = element_text(face = "bold"),
                   legend.box = "horizontal",
                   legend.spacing = unit(0,"mm"),
                   legend.box.background = element_rect(color="black", size=1))

Code

ggsave("Bathy_Un140.png", width =6 ,height = 6.5, units = "in")

Here I cut a shapefile to the 140m depth. This was the original depth I was going with before exploring the species depths but can customize once we agree on a depth.

Code

bath_cut<-bath_fort[bath_fort$z<0 & bath_fort$z> -140,]

#Rasterize
bath_sp<-SpatialPoints(coords = cbind(bath_cut$x, bath_cut$y))
r<-raster(ncols=297, nrows=297, xmn = -82,xmx = -62, ymn = 24.9, ymx = 45)
bath_r<-rasterize(bath_sp,r)
polys<-rasterToPolygons(bath_r)
poly_sf<-st_as_sf(polys)
poly_sf$id<-c(1:nrow(poly_sf))

#Remove grid points to isolate unwanted area
toremove<-c(10,12,14,16,9,11,13,39,41,43,40,42,44)
poly_sf<-poly_sf[!poly_sf$id %in% toremove,]

#Get Multiple Polygons to subset to 1 needed
poly_un<-st_union(poly_sf)
poly_cast<-st_cast(poly_un, "POLYGON")
poly_sf<-st_as_sf(poly_cast)
poly_fort<-fortify(poly_sf)
poly_fort$id<-c(1:nrow(poly_fort))
poly_fort<-poly_fort[poly_fort$id==4,]

#Buffer the final poly
poly<-poly_fort
# poly<-st_simplify(poly,dTolerance = .05)
# poly<-st_buffer(poly, dist = 0.1, nQuadSegs = 40)

ggplot(poly)+
  geom_sf(fill="#3788A6",color="#281E37", linewidth=0.75)+
  # geom_sf(data=poly_fort,  color="black", fill="black",alpha=0.3)+
  geom_sf(data=land_atl, inherit.aes = F, linewidth=0.5)+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  theme_minimal()+theme(panel.border = element_rect(color="black",fill=NA))

Code

ggsave("PolyTest.png", height = 9, width = 7, units = "in")

I cut the bathymetric data to the under 140m depth so that we could get a better view of the depth within our possible new area of analysis. As well as showing the orginal bathymetric map with the proposed area of analysis at 140m.

Code

#crop & mask
bath_rast<-as.raster(bathy)
bath_crop<-crop(bath_rast, extent(poly))
bath_mask<-mask(bath_crop, poly)

# bath_fort<-fortify(bath_mask)
bath_df<-as.data.frame(bath_mask, xy=T)
bath_df<-na.omit(bath_df)

# bath_df<-bath_df[bath_df$layer> -150,]

ggplot(bath_df)+
  geom_raster(aes(x=x, y=y, fill=layer))+
  geom_sf(data=land_atl, inherit.aes = F, linewidth=0.5)+
  geom_sf(data=poly, color="black", fill=NA, linewidth=0.75)+
  scale_fill_viridis_c(option = "G")+
  labs(title="Bathymetry ~ < 140m", x=NULL, y=NULL, fill="Depth (m)")+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  coord_sf(xlim = c(-82,-65), ylim = c(24.9,45))+
  theme(panel.border = element_rect(color="black",fill=NA),
        panel.background = element_rect(fill="white"),
        legend.position = c(0.83,0.17),
        title = element_text(face="bold", size = 13))

Code

ggsave("Bathy_CutShape.png", width = 5, height = 6.5, units = "in")

#All Bathy w/ shape outline
ggplot(bath_fort, aes(x=x, y=y, z=z))+
  geom_raster(aes(fill=z))+
  geom_sf(data=poly, inherit.aes = F, color="black", fill=NA, linewidth=0.75)+
  geom_sf(data=land_atl, inherit.aes = F, fill="gray70")+
  scale_fill_viridis_c(option = "G")+
  labs(title="Bathymetry < 140m", x=NULL, y=NULL, fill="Depth(m)")+
  scale_x_continuous(expand = c(0,0))+
  scale_y_continuous(expand = c(0,0))+
  coord_sf(xlim = c(-82,-62), ylim = c(24.9,45))+
  theme_bw()+theme(legend.position = c(0.912,0.135),
                   title = element_text(face = "bold"))

Code

ggsave("Bathy_wShape.png", width =6 ,height = 6.5, units = "in")